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Abstract. Several years of neutral particle measurements by the NASA/IBEX mis- 
sion have yielded direct observations of interstellar neutral helium and oxygen. The 
data indicate the presence of secondary neutral helium and oxygen, which are created 
within the heliosphere by charge exchange involving helium or oxygen ions. This con- 
tribution describes a detailed conserving calculation method based on Keplerian orbits 
that has been developed to characterize helium distribution functions throughout the 
heliosphere, in particular in the innermost heliosphere, while accounting for loss and 
production of neutral particles along their path. Coupled with global heliosphere mod- 
els of plasma distributions, this code is useful for predicting the fluxes of heavy neutral 
atoms at spacecraft detectors, so enabling inferences on the characteristics of the inter- 
stellar medium. 



1. Introduction 

On its galactic journey, the Sun traverses a cloud of interstellar gas that is moderately 
warm and dense, and is partially ionized, meaning that it also contains neutral atoms 
(dominantly hydrogen, but also helium, oxygen, and others). While the interstellar 
plasma interacts with the solar wind to form the global heliosphere, interstellar neutrals 
are unaffected and proceed into the heliosphere in their original velocity state. As such 
they are direct messengers of the inte rstellar medium, detectable by spacecraf t in the in- 
ner solar system (like NASA/IBEX; Bzo wski et alJl2012l : lM5bius et al.ll2012[) . This has 



intensified interest in accurate heliospheric neutral modeling, including secondaries. 
The entrance of hea vy interstellar neutrals into the heliosphere has been modeled 



in de tail in the past (e.g., iMiiller & ZarJdl2004l : llzmodenov et al.ll2004l : iBzowski et al. 



20121. and references therein). These calculations are typically variants of models and 
calculations of neutral interstellar hydrogen in the heliosphere, which have a long his- 
tory spanning four decades (for reviews, see lAxfordl 119721: IZardd 1 19991 : llzmodenov! 



l200rj) . Simulations treat the heavies as test particles on the background of the global hy- 
drogen/proton heliosphere. Monte-Carlo methods are useful in the global heliosphere, 
but statistics are poorer in the innermost heliosphere where IBEX measures. Reverse- 
trajectory methods are therefore popular to calculate heavy neutrals in the heliosphere, 
but typically the input flux is approximated at an outer boundary that is still well inside 
the heliosphere. This paper focuses on a method combining the best of both approaches, 
starting from the pristine interstellar medium, but calculating the phase space distribu- 
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tion function (PDF) at an arbitrary point of interest in a conserved way through reverse 
trajectories while allowing for arbitrary ionization (loss) and production terms. 



2. Scientific framework and numerical details 

In a heliocentric frame of reference, the motion of the neutrals is governed by the central 
gravitational potential of the Sun, such that the neutrals are on Keplerian trajectories, 
their paths deflected close to the Sun. The potential is weak for hydrogen which ex- 
periences an additional outward radiation pressure force that balances most of gravity. 
During solar maximum when radiation levels are highest in the solar cycle, radiation 
even overcompensates gravity, leading to an effectively repulsive potential for hydrogen 
during these times. Heavy neutrals are easier to handle, and the unimportance of radi- 
ation pressure for them means that the corresponding central potential is conservative 
and time-independent. 

Primary (interstellar) particles are on open trajectories (hyperbolae, and possibly 
parabolae as limiting case). Each individual trajectory can be described geometrically 
by specifying the orbital plane, the direction of perihelion, and the orbital eccentricity 
e. If 9, the angle that a particle's position vector makes with the direction of perihelion, 
is given, then the particle's radial distance and velocity are determined with simple 
formulae from celestial mechanics (r(8); v r {9); Vg{&)); the time to or since perihelion 
t{0) can also be obtained relatively easy albeit not in closed form. If any of the latter 
quantities is given instead, the other variables are determined just as easy from it, even if 
sometimes two solutions arise describing situations which are symmetric with respect to 
the perihelion (e.g., 6{r)). Each individual trajectory can also alternatively be described 
physically by specifying conserved quantities including total specific energy E, specific 
angular momentum /, and the eccentricity vector a dMuller & Coher]|2012l) . Again with 
trajectories defined in this manner, any one given value for angle, distance, velocity 
or time will yield all the respective others. There are of course numerous relations 
that link orbital elements to the physical quantities (for instance, the orbital plane is 
perpendicular to /, and e and the perihelion distance r m! „ depend only on / and E). 

On their way through the heliosphere, atoms have a probability of being lost. In 
the case of interstellar helium, in order of importance for helium atoms detected at a 1 
AU distance (e.g., location of IBEX), the processes are loss due to photoionization, loss 
due to charge exchange (ex.) with slow helium ions in the outer helios heath, double 



ex. w ith solar wind alpha particles, and ex. with solar wind protons (Bz owski et al. 
20121 ). Earlier studies (e.g., iMuller & Zankl f2004) were limited in admitting only ex. 
with protons as a loss process. Secondary helium is created as the neutral product 
of a ex. collision between a helium ion and a neutral, the dominant interaction being 
bow-shock decelerated interstellar He + exchanging charge with interstellar helium in 
the outer heliosheath (loss of a primary particle and at the same time creation of a sec- 
ondary neutral in a different velocity state, reflecting the underlying plasma distribution 
function). Of course, once created, a secondary neutral on its further trajectory suffers 
the same types of losses as the primary neutrals. 

Two frames of reference are most appropriate for the neutral calculations: (1) The 
"ISM frame" in which the x axis points antiparallel to the interstellar (ISM ) flow so tha t 
the ISM flow vector at infinity is (u xJSM , 0, 0) with u x j S m - -26.3 km/s dWittell2004h : 
and (2) the "Kepler frame" in which the x axis points to the perihelion of the trajectory 
in question, and the x-y plane is the orbital plane. For each trajectory, a z -l x = L = in 
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a)X=-10, Y=0 C)X=-3.5, Y=2 
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Figure 1. Phase space densities of primary neutral interstellar helium at down- 
wind locations, on the symmetry axis ((-10.0, 0.0), panels (a) and (b)), and at a 
location slightly off the symmetry axis at (-3.5, 2.0), distance 4.0 AU (panels (c) and 
(d)). Both locations are in the helium focusing cone, with the PSD being a deformed 
torus loosely perpendicular to the x-axis; see text for details. Presented for each PSD 
are two perpendicular cuts through the 3D object, with the cut plane specified (panel 
(d) and its inset have different cut planes, as indicated). The maximum PSD value in 
the ISM is normalized to 1 (logarithmic contour levels with 3 lines per dex. 



its Kepler frame, leading to a convenient reduction of conserved quantities that need to 
be considered. Rotation transformations relate ISM frame coordinates to Kepler frame 
coordinates, and also to other coordinates of choice such as a heliocentric frame of 
reference associated with ecliptic longitude and latitude. Spherical polar coordinates 
like the latter are typically used for the heliosphere plasma background model, r is 
invariant under these transformations. 

The strategy to calculate a phase space distribution function (PDF) at a point of 
interest ro is essentially to investigate all possible trajectories that pass through ro. For 
the primary PDF some of these trajectories have to be disregarded because they are 
not connected to the ISM (everything outside of the heliospheric bow shock); similarly, 
some trajectories that do not intersect with the outer heliosheath do not contribute to- 
wards the secondary PDF. In principle, all these trajectories passing through r*o can be 
obtained by scanning the entire 3D velocity space; for each individual phase space lo- 
cation (ro, vq), the reverse trajectory is uniquely determined and the velocity at infinity, 
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or at any point where the trajectory intersects the outer heliosheath, is calculated in one 
algebraic step, while preserving the conserved quantities by design. 

Neglecting secondary neutrals for the moment, the prescription for the primary 
PDF is then as follows: By Liouville's theorem, the raw phase space density at (ro, vo) 
is the same as the value of the ISM Maxwellian at the calculated velocity at infinity. 
Attenuation factors arise due to photoionization and charge exchange; the photoioniza- 
tion factor can be cast into an analytic form involving essentially only expressions for 
9 and for the asymptotic #00 at infinity, the latter a function of e. As we want to take 
the plasma background from an external global heliosphere model, the ex. attenuation 
requires a ex. loss integration along the trajectory from infinity (in reality, from the 
bow shock crossing) to ro, which is carried out at the r-grid nodes of the background 
model, with 6{r) determined algebraically (conserving by design) and the transforma- 
tion matrix between Kepler and global heliosphere frames being constant during this 
calculation. 

In practice, instead of 3D scanning, the locations of local maxima in the PDF are 
predicted algebraically: There are typically two maxima associated with the primary 
neutrals, called direct- and indirect-path vel ocities, where the indi rect path is longer 
and leads the particle closer to the Sun (e.g., iMuller & Cohenll2012l) . These velocities 
can be obtained in one step in the ISM frame by requiring the velocity at infinity to be 
UxiSM, which is the location of the maximum of the ISM Maxwellian: 



fu y ~ u x ism 
(v x , v y ) = (u xJ sm, 0) -(y, r - x) ; l z = -r 

W z 2 



1± Jl+4. U r °~ X 



u Iism y 2 



(i) 



with (x,y) - ?o, f M the constant of gravity multiplied by the solar mass, and the upper 
sign the direct path. The derivation typ ically makes u se of the conserved quantities; 
similar results have been summarized by lAxfordl (119721) and others. 

The simplest practical strategy for the primary helium PDF is hence to obtain the 
location of the two maxima with equation (1), and then scan velocity space around them 
until the phase space density has dropped off by, say, 4 orders of magnitude which is 
equivalent to ~4 thermal velocities, beyond which it can safely be neglected. A more 
sophisticated procedure is to develop the equivalent of eq. (1) for the fixed ISM frame 
with arbitrary interstellar vectors it and with it, map the locus of the edge of the ISM 
Maxwellian to velocities at ro. This gives a closed velocity surface at ro within which 
the non-negligible PDF resides. 



3. Primary helium results 

Figure 1 sho ws two examples of applying the code, us ing typical I SM characteristics 
( Wittej l2004h . a typical heliospheric background model (Muller et al.ll2006h . and typical 



photoionization and ex. rates. In the helium focusing cone on the downwind symmetry 
axis (Fig. la, b), the primary helium PSD is degenerate, with all planes containing 
the symmetry axis being equally preferred orbital planes for the interstellar particles. 
This leads to a torus-like PDF, with the torus perpendicular to the v x axis and the cross 
section of the torus slightly gravitationally deformed away from a circular shape. The 
locus of maxima (64% of the interstellar Maxwellian maximum due to losses) is a circle 
with radius J-2f M /x = 13.3 km/s at a constant v x = u x jsm- 
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a) direct b) indirect c) indirect 




Figure 2. Phase space density cuts of primary neutral interstellar helium at a 
sidewind location (0, 2) through the direct (a) and the indirect (b, c) PDF. 



Moving fo slightly off the symmetry axis (Fig. lc, d), it can be seen that the degen- 
eracy is broken. The plane of the torus is now tilted with respect to the v y -v z plane, and 
the torus now connects a thick knot (maximum 62%) to a small knot (34%), with the 
connections having only very low phase space density. Going further off-axis makes 
it clear that the thick knot is associated with direct path trajectories; the small knot is 
associated with the indirect path (more prone to losses and hence less dense), and they 
are no longer connected. Inspecting upwind directions, the direct portion approaches a 
Maxwellian, whereas the indirect portion is severely deformed, having to travel to the 
Sun first and then back to the point of interest, with the incumbent losses along the way. 

An intermediate situation is shown in Figure 2, for the sidewind location (0, 2). 
Direct and indirect PDF are well separated. The direct peak (64%) occurs for negative 
v x and negative v y , whereas both are positive for the indirect peak (5%). The 3D shapes 
are lenticular for the direct PDF, and even more irregular for the indirect one. The 
velocity volume occupied by these PDF traces back to the ISM and its Maxwellian. 
In configuration space, all the associated trajectories go through the location (0, 2), 
but at a remote reference distance, for example a sphere of radius 1000 AU, they fill a 
certain area. Note that not all interstellar particles passing through that area intersect 
the location (0, 2); in that sense, one can only talk about mapping the ISM Maxwellian 
to the PDF at (0, 2) in velocity space, but not in configuration space. 

With the fast, accurate characterization of the neutral PDF in hand, it is straight- 
forward to perform integrals in velocity space including moments of the PDF, to arrive 
at configuration space maps of number density, bulk velocity, temperature, loss and 
production rate, and others. All of these quantities also can be obtained separately for 
direct and indirect neutrals. Given the secondary production rate, through forward cal- 
culation the peak of the secondary PDF at an arbitrary location fo can be constrained, 
making the calculation of the secondary PDF (scanning a velocity neighborhood of the 
expected center) more efficient. As mentioned above, for each individual secondary 
PDF point, an integration (nominally from bow shock to heliopause) is necessary to 
sum up the production terms and convolve the integration as usual with losses on the 
way to ro. This integration involves the knowledge of the primary PDF, and as such, the 
process has a recursive flavor to it. As seen from the vantage of heliospheric physics, 
one of the less resolved ingredients in the problem of secondary PDF is how the he- 
lium ions necessary for ex. are being handled correctly. This problem is outside the 
purview of the neutral solver exhibited in this contribution, but pertains more to the 
correct handling of heavy ions through the global plasma heliosphere models. 
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4. Conclusions 

Using the Kepler equations of celestial mechanics is an elegant, fast way to calculate 
trajectories of neutral heavy atoms originating in the local interstellar medium. The 
associated conserved quantities enable one-step calculations of particle locations and 
velocities that are accurate by design. The presented examples for helium show the 
sometimes complicated heliospheric phase space distribution functions. They empha- 
size the role of gravitational deflection and the effect of loss and production processes 
on the helium phase space distribution. The primary neutral helium PDF can be calcu- 
lated in this way at any arbitrary point in the heliosphere, and its various moments yield 
maps of effective temperature, bulk velocity, and density. The latter are elevated in the 
region of the helium focusing cone, which is also where direct and indirect solutions 
merge together and therefore become equally important. 

The developed computational met hod, if app lied only to primaries, is somewhat 
parallel to the hot model of neutrals dFahrl [l97lT) . Yet the same method applies to 
secondary helium neutrals as well, which is a novel feature. The method furthermore 
treats the plasma as background and therefore can be used with a variety of global helio- 
spheric models to study the sensitivity of the helium signal to heliospheric asymmetries 
or other differences in the global models. Lastly, the mathematics of the method are not 
affected by a time-dependence neither of the background plasma nor of the photoion- 
ization rate; the method is therefore adaptable to realistic, time-dependent scenarios, 
necessitating only an increased level of housekeeping to match times in the trajectory 
with the time-dependent inputs. 
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